function [C_f f]=Ford_Fullerson(C);
 %[C_f f]=Ford_Fullerson(C)
 %C:capacity of the network
 %C_f:residual network
 %f:flow
 
C_size=size(C);
f=zeros(C_size);
C_f=C;
%
while(1)
 %get a path
 P=find_path(C_f);%P=[u,v;]
if isempty(P)
    break;
end
c_f=min(C_f(P(:,1)+(P(:,2)*C_size(2)-C_size(2))));%flow of P
f_c=zeros(C_size);
f_c(P(:,1)+(P(:,2)*C_size(2)-C_size(2)))=c_f;%cloumn
 %f=f+f'
 f=f+f_c;
 %residual graph
 C_f(P(:,1)+(P(:,2)*C_size(2)-C_size(2)))=C_f(P(:,1)+(P(:,2)*C_size(2)-C_size(2)))-c_f;
 C_f(P(:,2)+(P(:,1)*C_size(2)-C_size(2)))=C_f(P(:,2)+(P(:,1)*C_size(2)-C_size(2)))+c_f;
end
end

function P=find_path(C);
C_size=size(C);
if C_size(1)<2 || C_size(1)~=C_size(2)
    P=0;
    return;
end
L=C_size(1);
S=1;
T=2:L;
P_record=[];
while 1
    [max_value1 max_idx1]=max(C(S,T),[],1);
    [max_value2 max_idx2]=max(max_value1,[],2);
    max_value=max_value2;
    max_idx=[max_idx1(max_idx2) max_idx2];
    max_path=[S(max_idx(1)) T(max_idx(2))];
    if max_value==0
        P=[];
        return;
    end
    %[min_value min_idx]=min(C.* (C(S,T)>0));
    P_record=[P_record;max_path];
    if max_path(2)==L
        break;
    end
    S=[S max_path(2)];
    T=[T(1:max_idx(2)-1),T(max_idx(2)+1:length(T))];
end
%back for path
P=[];
target=L;
for idx=size(P_record,1):-1:1
    if P_record(idx,2)==target
        P=[P_record(idx,:); P];
        target=P_record(idx,1);
    end
    if target==1
        break;
    end
end
end

